/* * * * M A P V I S I O N   L I C E N S E  *****************************************************************
 *
 * This MapVision software is © Filip Zavadil 2013. It is licensed under the LGPL license.
 * This MapVision software dynamically links to unmodified Qt5 Library.
 *
 * The Qt5 Library is © 2013 Qt® and/or its subsidiary(-ies), and is licensed under the
 * GNU Lesser General Public License version 2.1 with Qt® LGPL exception version 1.1.
 * Qt5 library is free software; you can redistribute it and/or modify it under the terms of the
 * GNU Lesser General Public License, version 2.1, as published by the Free Software Foundation.
 * Qt5 library is provided “AS IS”, without WARRANTIES OR CONDITIONS OF ANY KIND, EITHER EXPRESS OR
 * IMPLIED INCLUDING, WITHOUT LIMITATION, ANY WARRANTIES OR CONDITIONS OF TITLE, NON-INFRINGEMENT,
 * MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
 *
 * As an additional permission to the GNU Lesser General Public License version 2.1, the object
 * code form of a “work that uses the Library” may incorporate material from a header file that is
 * part of the Library. You may distribute such object code under terms of your choice,
 * provided that: (i) the header files of the Library have not been modified; and (ii)
 * the incorporated material is limited to numerical parameters, data structure layouts, accessors,
 * macros, inline functions and templates; and (iii) you comply with the terms of Section 6 of the
 * GNU Lesser General Public License version 2.1.
 *
 * Moreover, you may apply this exception to a modified version of the Library, provided that such
 * modification does not involve copying material from the Library into the modified Library’s header
 * files unless such material is limited to (i) numerical parameters; (ii) data structure layouts;
 * (iii) accessors; and (iv) small macros, templates and inline functions of five lines or less in length.
 *
 * Furthermore, you are not required to apply this additional permission to a modified version of the Library.
 *
 * You should have received a copy of the GNU Lesser General Public License along with this package;
 * if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA

 * The source code for Qt 5.X SDK is available from Qt Project here:
 * http://download.qt-project.org/official_releases/qtcreator/2.7/2.7.0/qt-creator-2.7.0-src.zip.
 *
 * ****************************** M A P V I S I O N    S O F T W A R E  ************************************
 *
 * Author:       Filip Zavadil
 * Description:  Software for displaying and exporting map tiles.
 * The sources and executables of MapVision software are available:
 *    - from project hosting on CodePlex.com :   http://mapvision.codeplex.com
 *    - from GIT repositories on Gitorious.org:  https://gitorious.org/mapvision
 * Visit MapVision fanpage on Facebook.com:      https://www.facebook.com/pages/Mapvision/631827753513935
 *
 * *********************************************************************************************************
*/


#include "zavaproj.h"
#include <QString>
#include <QStringList>
#include <QDebug>
#include <QTimer>
#include <QEventLoop>
#include <QMessageBox>
#include <QWidget>
#include "proj_api.h"

ZavaProj::ZavaProj()
{

}

QPoint ZavaProj::lonlat2tileFormat(QPointF lonlat,int zoom,Yaxis yaxis)
{
    QPoint xy;
    double lon,lat;
    lon=lonlat.x();    lat=lonlat.y();

    int m=pow(2.0, zoom);
    xy.setX((int)(floor((lon + 180.0) / 360.0 * m)));
    xy.setY((int)(floor((1.0 - log( tan(lat * M_PI/180.0) + 1.0 / cos(lat * M_PI/180.0)) / M_PI) / 2.0 * m)));
    // conversion to TMS
    if (yaxis==BOTTOM) xy.setY(m-xy.y()-1);
    return xy;
}


QPoint ZavaProj::lonlat2tileFormat(double lon, double lat,int zoom,Yaxis yaxis)
{
    QPoint xy;
    int m=pow(2.0, zoom);
    xy.setX((int)(floor((lon + 180.0) / 360.0 * m)));
    xy.setY((int)(floor((1.0 - log( tan(lat * M_PI/180.0) + 1.0
                                    / cos(lat * M_PI/180.0)) / M_PI) / 2.0 * m)));
    // conversion to TMS
    if (yaxis==BOTTOM) xy.setY(m-xy.y()-1);
    return xy;
}

QPointF ZavaProj::tileFormat2lonlat(QPoint xy, int zoom, Yaxis yaxis)
{
    QPointF lonlat;
    double x,y;

    x=xy.x();    y=xy.y();
    // conversion to TMS
    int m=pow(2.0, zoom);
    if (yaxis==BOTTOM) y = m-y-1;
    lonlat.setX(x / m * 360.0 - 180.0);
    double n = M_PI - 2.0 * M_PI * y / m;
    lonlat.setY( 180.0 / M_PI * atan(0.5 * (exp(n) - exp(-n))));
    return lonlat;
}

QPointF ZavaProj::tileFormat2lonlat(int x, int y,int zoom,Yaxis yaxis)
{
    QPointF lonlat;
    // conversion to TMS
    int m=pow(2.0, zoom);
    if (yaxis==BOTTOM) y = m-y-1;
    lonlat.setX(x / m * 360.0 - 180.0);
    double n = M_PI - 2.0 * M_PI * y / m;
    lonlat.setY( 180.0 / M_PI * atan(0.5 * (exp(n) - exp(-n))));
    return lonlat;
}

QPointF ZavaProj::pixels2lonlat(QPoint pixelPos, int zoom)
{
    double pixelX=pixelPos.x();
    double pixelY=pixelPos.y();
    double mapSize = 256 << zoom;
    double x = (Clip(double(pixelX), 0, mapSize - 1) / double(mapSize)) - 0.5;
    double y = 0.5 - (Clip(double(pixelY), 0, mapSize - 1) / double(mapSize));

    double lat = 90 - double(360) * atan(exp(-y * double(2) * M_PI)) / double(M_PI);
    double lon = double(360.0) * x;

    return QPointF(lon,lat);
}

QPoint ZavaProj::lonlat2pixels (QPointF lonlat, int zoom)
{
    QPointF min_maxLon;QPointF min_maxLat;
    double lon=lonlat.x();
    double lat=lonlat.y();
    double minLat=min_maxLat.x();
    double maxLat=min_maxLat.y();
    double minLon=min_maxLon.x();
    double maxLon=min_maxLon.y();


    minLat=-85.05112878;
    maxLat=85.05112878;
    minLon=-180;
    maxLon=180;



    lat = Clip(lat, minLat, maxLat);
    lon= Clip(lon, minLon, maxLon);

    double x = (lon + 180) / 360;
    double sinlat = sin(lat * M_PI / 180);
    double y = 0.5 - log((1 + sinlat) / (1 - sinlat)) / (4 * M_PI);

    uint mapSize = 256 << zoom;
    int pixelX = (int) Clip(x * mapSize + 0.5, 0, mapSize - 1);
    int pixelY = (int) Clip(y * mapSize + 0.5, 0, mapSize - 1);

    return QPoint(pixelX,pixelY);
}

QPoint ZavaProj::pixels2tileFormat(QPoint xy, int zoom, ZavaProj::Yaxis origin)
{
    QPointF lonlat=this->pixels2lonlat(xy,zoom);
    //  if (lonlat.x()>85.0511) lonlat.setX(85.0510);
    //   if (lonlat.y()>85.0511) lonlat.setY(85.0510);


    return this->lonlat2tileFormat(lonlat,zoom,origin);
}

QPoint ZavaProj::tileFormat2pixels(QPoint xy, int zoom, ZavaProj::Yaxis origin)
{
    QPointF lonlat = this->tileFormat2lonlat(xy,zoom,origin);
    if (lonlat.x()>85.0511) lonlat.setX(85.0510);
    if (lonlat.y()>85.0511) lonlat.setY(85.0510);
    return this->lonlat2pixels(lonlat,zoom);
}

int ZavaProj::convertYaxis(int yTile, int zoom)
{
    return pow(2,zoom)-yTile-1;
}

QPoint ZavaProj::convertYaxis(QPoint tile, int zoom)
{
    return QPoint(tile.x(),pow(2,zoom)-tile.y()-1);
}

QRect ZavaProj::rectForTile(QPoint tile)
{
    // returns TMS Y-origin (bottom)
    return QRect(tile.x()*256,tile.y()*256,256,256);

}

QPoint ZavaProj::pointForTile(QPoint tileFormat,int zoom)
{
    return QPoint(tileFormat.x()*256,this->convertYaxis(tileFormat.y(),zoom)*256);
}



double ZavaProj::Clip(double n, double minValue, double maxValue)
{
    return std::min(std::max(n, minValue), maxValue);
    //return qMin(qMax(n, minValue), maxValue);
}

double ZavaProj::dms2double(QString input){

    bool ok=false;

    bool negative=false;

    double output,deg,min,sec=0;

    // remove all white spaces
    while (input.contains(" "))    {
        input=input.remove(" ");
    }

    // different formating
    input=input.replace(",",".");
    input=input.replace("°","°");
    input=input.replace("'","'");
    input=input.replace("”","\"");

    if (input.contains("S") || input.contains("W")){
        negative=true;
    }
    input=input.remove("S");
    input=input.remove("N");
    input=input.remove("E");
    input=input.remove("W");


    // DDD°MM'SS.SS" - sexagecimal
    if (input.contains("°") && input.contains("'") && input.contains("\"")){

        deg= input.mid(0,input.indexOf("°")).toDouble(&ok);
        if (ok) min=input.mid(input.indexOf("°")+1,input.indexOf("'")-input.indexOf("°")-1).toDouble(&ok);
        if (ok) sec=input.mid(input.indexOf("'")+1,input.indexOf("\"")-input.indexOf("'")-1).toDouble(&ok);
        output=deg;
        output+=(min/60);
        output+=(sec/3600);

        if (negative) output=-output;
        if (!ok) output = NAN;
        return output;
    }
    else if (input.contains("°") && input.contains("'") ){
        // DDD°MM.MMM' - sexagesimal
        deg= input.mid(0,input.indexOf("°")).toDouble(&ok);
        if (ok) min=input.mid(input.indexOf("°")+1,input.indexOf("'")-input.indexOf("°")-1).toDouble(&ok);
        output=deg;
        output+=(min/60);

        if (negative) output=-output;
        if (!ok) output = NAN;
        return output;
    }
    else if (input.contains("°")){
        // DDD.DDDDD° - sexagesimal
        input=input.remove("°");
        output=input.toDouble(&ok);
        if (negative) output=-output;
        if (!ok) output = NAN;
        return output;
    }
    else
    {
        // decimal format
        output=input.toDouble(&ok);
        if (negative) output=-output;
        if (!ok) output = NAN;
        return output;
    }

}

QPointF ZavaProj::dms2lonlat(QString input)
{

    if (input.count(",")==1){
        QStringList coors=input.split(",");
        if (coors[0].contains("W")||coors[0].contains("E"))
        {
            return QPointF(this->dms2double(coors[0]),this->dms2double(coors[1]));
        }
        else
        {
            return QPointF(this->dms2double(coors[1]),this->dms2double(coors[0]));
        }
    }
    else if (input.count(" ")==1){
        QStringList coors=input.split(" ");
        if (coors[0].contains("W")||coors[0].contains("E"))
        {
            return QPointF(this->dms2double(coors[0]),this->dms2double(coors[1]));
        }
        else
        {
            return QPointF(this->dms2double(coors[1]),this->dms2double(coors[0]));
        }
    }
    else
    {
        return QPointF(NAN,NAN);
    }
}


QString ZavaProj::double2dms(double input, int prec)
{

    double positive_degrees,degrees,decimal_part,minutes,seconds,remains= 0;

    positive_degrees=std::abs(input);

    degrees=int(positive_degrees);
    decimal_part=60*(positive_degrees-degrees);
    minutes=int(decimal_part);
    decimal_part=60*(decimal_part-minutes);
    seconds=int(decimal_part);
    remains=(decimal_part-seconds);

    return QString::number(degrees,'f',0)+"°"+
            QString::number(minutes,'f',0)+"'"+
            QString::number(seconds+remains,'f',prec)+"\"";

}

QString ZavaProj::lonlat2dms(QPointF lonLat, int prec)
{
    QString lon,lat,DMS;
    lon=this->double2dms(lonLat.x(),prec);
    lat=this->double2dms(lonLat.y(),prec);
    DMS = (lonLat.y()>0) ? lat + "N" : lat + "S";
    DMS += ", ";
    DMS +=(lonLat.x()>0) ? lon + "E" : lon + "W";
    return DMS;
}

double ZavaProj::getDistance(QPointF lonlatA, QPointF lonlatB)
{
    //http://www.movable-type.co.uk/scripts/latlong.html
    int R = 6371;
    double dLat = lonlatA.y() - lonlatB.y();
    double dLon = lonlatA.x() - lonlatB.x();
    double rad =  M_PI / double(180);
    dLat *= rad;
    dLon *= rad;

    double a = sin(dLat/2) * sin(dLat/2) + sin(dLon/2) * sin(dLon/2) * cos(lonlatA.y()*rad) * cos(lonlatB.y()*rad);
    double c = 2 * atan2(sqrt(a),sqrt(1-a));

    return R * c ;
}

QPointF ZavaProj::project(QPointF lonlat, QString proj4CODE, bool isDegree)
{
    double x,y;
    projPJ out,in;
    x=lonlat.x();
    y=lonlat.y();

    x*=DEG_TO_RAD;
    y*=DEG_TO_RAD;

    if (!(out = pj_init_plus(proj4CODE.toUtf8()))){
        QMessageBox msgBox;
        msgBox.setText("The input projection is invalid! ");
        msgBox.exec();
        return QPointF(NULL,NULL);

    }
    if (!(in = pj_init_plus("+proj=latlong +ellps=WGS84 +datum=WGS84 +no_defs"))) return QPointF(NULL,NULL);


    if (!pj_transform(in,out,1,1,&x,&y,NULL)) {
        if (isDegree){
            x*=RAD_TO_DEG;
            y*=RAD_TO_DEG;
        }
        return QPointF(x,y);
    }
    else return QPointF(NULL,NULL);
}

QPointF ZavaProj::lonlat2UTM(QPointF lonlat)
{
    double x = lonlat.x() * 20037508.34 / 180;
    double y = log(tan((90 + lonlat.y()) * M_PI / 360)) / (M_PI / 180);
    y = y * 20037508.34 / 180;
    return QPointF(x,y);
}

QPointF ZavaProj::UTM2lonlat(QPointF xy)
{
    double lon = (xy.x() / 20037508.34) * 180;
    double lat = (xy.y() / 20037508.34) * 180;

    lat = 180/M_PI * (2 * atan(exp(lat * M_PI / 180)) - M_PI / 2);
    return QPointF(lon, lat);
}

QString ZavaProj::lonlat2dms(QPointF lonLat)
{
    QString lon,lat,DMS;
    lon=this->double2dms(lonLat.x(),2);
    lat=this->double2dms(lonLat.y(),2);
    DMS = (lonLat.y()>0) ? lat + "N" : lat + "S";
    DMS += ", ";
    DMS +=(lonLat.x()>0) ? lon + "E" : lon + "W";
    return DMS;
}



